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1. Introduction 

During a floating zone experiment involving the growth of indium on a recent orbiter 
mission, (STS 32) oscillation of the zone shapes were observed to occur in response to the 
background acceleration. An understanding of the nature of the response of the zone shape to 
forced (g-jitter) oscillations and predictions of its impact on future experiments is of great 
interest not only to the Pi’s of but to other commercial and academic investigators who plan to 
fly similar experiments in the orbiter and on space station.Motivated by this, a 15 month study 
was undertaken to analyze the nature of the g-sensitivity of the STS 32 floating zone crystal 
growth experiment. Numerical models were used to describe the time-dependent free surface 
motion of the zone as it responds to the spacecraft residual acceleration. Relevant experimental 
data concerning the acceleration environment was obtained from the Honeywell in Space 
Accelerometer (HISA) investigators through MSFC's ACAP program and processed and 
analyzed. For the indium floating zone experiment a series of calculations were made using time- 
dependent axial accelerations g(t). The form of g(t) included simple sinusoidal disturbances as 
well as actual data (subject to appropriate filtering) measured on the STS 32 mission. The study 
focused the calculation of the response of the free surface of the zone as well as the internal 
flows and internal heat transfer. The influence of solidification on the response of the zone shape 
was also examined but found to be negligible. 

The results are summarized in section 3 and details are supplied in the appendix. In 
section 4, published papers and presentations given that involve this work are listed. References 
are given in section 5. 

2. Background 

Zhang and Alexander [ 1 ] have approached the problem of determining the axisymmetric 
response of the shape of the free surface of a cylindrical liquid column bounded by two solid 
regions is modelled by 1-D system of non-linear equations. It is found that the sensitivity of the 
zone shape depends on the static Bond number, B 0 , aspect ratio and viscosity as well as the 
amplitude and frequency of the disturbance. The general trend is an increase in tolerable residual 
gravity with increasing frequency. At the eigenfrequencies of the zone, however, there are 
dramatic deviations from this trend. At these frequencies the tolerable residual gravity level can 
be two orders of magnitude lower at this frequency. Aspect ratios ranging from 96.3 to 82.7% of 
the Rayleigh limit were examined. For these cases, the frequencies associated with the lowest 
tolerable acceleration have been found in the 10-2-1 O' 1 Hz range. In terms of previously recorded 
and predicted residual accelerations the sensitive frequency ranges for the cases examined are 10- 
2 -10- 1 Hz and 1-10 Hz. Maximum tolerable residual gravity levels as low as 10- 6 g have been 
calculated. The effect of viscosity is seen to increase the tolerable acceleration level for all 
frequencies. The equilibrium shape, as determined by the steady background acceleration, has a 
pronounced effect at low frequencies. A change in slenderness of the bridge markedly changes 
the sensitivity to residual acceleration as the Rayleigh limit is approached. 

To asses the influence of g-jitter on thermocapillary flow it is useful to consider the 
dimensionless groups involved. The dimensionless groups associated with mixed 
thermocapillary convection and buoyancy-driven convection are [2]: 

The Dynamic Bond number: 

The Marangoni number: 

The Grashof number 

The Prandtl number 

The Surface Reynolds number 


Bd = pgL 2 /y T 
Ma = y-jATLAik 
G r = ATpgL^/v 2 
Pr = v/k 
R s = MaPr 1 , 
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Here p is the liquid density, y T = dy/dT, p and v are the shear and kinematic viscosities of the 

fluid melt, g is the acceleration, L is some characteristic length scale, AT is a characteristic 
temperature difference and k is the liquid's thermal diffusivity. 

The relative importance of gravity and thermocapillary forces can be estimated as follows 

[2,3]: 

First, the velocity reference scale is obtained from a consideration of the balance of force 
tangent to the free surface [2] 





(i) 


If the velocity and temperature are scaled by V and AT, and lengths are scaled by L, then 
(1) can be recast as 


3u_ /mat] ae 
ax ^ dz ’ 


where u and 0 represent the dimensionless velocity and temperature, and x and z are the 
dimensionless lengths. 

For terms on both sides of (2) to be of the same order the magnitude of V* must be [2,3] 


y* _ NAT 
M- 


(3) 


The balance of linear momentum and energy in a Boussinesq fluid are given by 

a? i ~ _ ~ _ 

— + (gradv)v = -J-grad p + vAv + PjTg(t). 

at P 


(4) 


For steady flows Ostrach [2] used (3) to scale (4). For an aspect ratio of 1 and with R s «1 he 
obtained: 

p 

Av=-gradp + -^g(t), (5) 

where the pressure is scaled with |iV7L. 

For R s »1 a boundary layer scale, 8 = L/Vr s , was employed to scale the coordinate 
perpendicular to the free surface. For this steady boundary layer situation the component of (4) 
perpendicular to the interface has the (dimensionless) form 


dv dv 


= r 2/3^P + _1_^ | d 2y 

dy rP a x 2 


3y 2 + ir 8(,) ' 


( 6 ) 


At small values of R s the relative contribution of the buoyancy force to convection can be 
estimated from the ratio G r /R s • For large values of R s , the relative importance of buoyancy is 
found upon comparing the buoyancy term with the pressure gradient and is expressed via the 
ratio Gr/R s 2/3 . 

In order to estimate the effect of time-dependent periodic accelerations on a 
thermocapillary flow we adopt the approach of Langbein and Tiby and express the magnitude of 

the velocity caused by an acceleration occurring with a circular frequency co as 
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Strouhal Number (St) 


Fig. 1 Sensitivity regimes for thermocapillary flow expressed in terms of the Grashof (Gr), Strouhal (St) and surface 

Reynolds (Rs) numbers. After Alexander [3]. 

For thermocapillary flow to dominate, the ratio of Vg to V must be less than one. For 
this to occur the magnitude of the periodic acceleration must satisfy: 

V,. L ^-, (7) 

( 0 ) 2 + — )' /2 

L 4 

Equation (V) can be re-interpreted in terms of the ratio of the Grashof number, Gr, and the 
surface Reynolds number R s . The condition that buoyancy and thermocapillary forces are of the 
same magnitude is [3] 

G r = R s (S t 2 + 1) 1/2 , (8) 

where S,= toL 2 /v is the Strouhal number. This condition is illustrated graphically in Fig. 1. 
Three regimes of interest are identified. These regions are defined by the relative values of the 
Grashof, Strouhal and surface Reynolds numbers. In region I, the Strouhal number S t is less than 
one, and the condition that thermocapillary forces dominate is equivalent to that for a steady 
flow. In region II, S t is greater than one. In this region, for a fixed value of R s , the value of Gr 
required to give buoyancy forces equal weight to thermocapillary forces increases with 
increasing values of the Strouhal number. This essentially reflects the fact that the characteristic 
time for the fluid response greatly exceeds the period of the disturbance. Thus, thermocapillary 
forces dominate at higher values of than in region I. In region III buoyancy forces predominate 

over thermocapillary forces. . 

Alexander [3] estimated the sensitivity of surface-tension-driven flow to residual 
acceleration for a number of different thermocapillary flow systems (under the assumption that 
the free surface would not deform). In general, these estimates suggest that periodic vibrations 
(with amplitudes and frequencies characteristic of spacelab) will generally not lead to significant 
buoyancy effects except in extreme cases. These estimates are confirmed by our numerical 
modelling. For cases where one ore more of the boundaries are deformable, the density 
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difference across that boundary determines the extent to which g-jitter induced motion of the 
surface will influence the internal flow in the liquid. For the case under investigation here, the 
interface motion can, as we shall see in the next section and in appendix I, dominate the flow 
(depending on the frequency and amplitude of the g-jitter). 

More recent work by Zhang and Alexander [4] involving full axisymmetric calculations 
indicates that the isothermal ID model give a reliable indication of the frequency dependence of 
the shape sensitivity for isothermal and some non-isothermal flows. The full numerical models 
are, however, essential for a description of the internal flow and heat transfer (especially in the 
non-isothermal cases). 

During a floating zone experiment involving the growth of indium on a recent orbiter 
mission STS 32, oscillation of the zone shapes were observed to occur in response to the 
background acceleration. In addition to these observed oscillations, zone lengths in excess of the 
often quoted Rayleigh limit [5,6] were recorded. The Rayleigh limit occurs when the length of 

the zone exceeds 2k times its radius and is associated with the tendency of the surface tension of 
the liquid to cause the zone to break up. Zone lengths in excess of this limit can be explained in 
several ways, including melt volumes in excess of the right circular cylindrical volume 
corresponding to the zone length and feed rod radius, by surface contamination resulting in zero 
surface tension, highly non-planar boundaries separating the solid indium form the indium melt 
(causing stabilization due to the large deviation of the zone shape from a circular cylinder), or by 
the stabilizing action of surface-driven (thermocapillary) flow [7]. Regardless of the mechanism 
by which these excessively long zones were achieved, the response of the zone shape to forced 
(g-jitter) oscillations is of great interest not only to the Pi’s of but to other commercial and 
academic investigators who plan to fly similar experiments in the orbiter and on space station. 


3. The model floating zone 

Governing equations 

A full numerical model of the axisymmetric free surface response of a non-isothermal liquid 
subject to axial vibration was employed. We assume that the surface tension of the melt varies 
linearly with temperature. Thus, the presence of a temperature gradient along the free surface of 
the liquid bridge gives rise to a thermocapillary flow. 

Consider a cylindrical liquid zone (see Fig. 2) held between two parallel coaxial circular 
rigid disks (radius = R 0 ) separated by a distance L. The liquid is a non-isothermal incompressible 
Newtonian fluid. The bridge is held between the disks by surface tension. The free surface of 
the bridge is a gas-liquid interface and is described by r=R(z,t). Each disk is maintained at a 
constant temperature T 0 . Surface heating is provided through a function, Too(z), which represents 
the heater profile. The heat transfer coefficient at the free surface is denoted by h. In addition, 
we make the assumptions that the residual acceleration is parallel to the cylinder axis and the 
velocity and temperature field and the deformation of free surface are axisymmetric. 
Furthermore, we let the surface tension at the free surface vary linearly with temperature and 
assume that the Boussinesq approximation holds. 

The governing equations are made dimensionless by scaling length, time and velocity 
with Ro, R 0 * and U* respectively. Here U* is a characteristic velocity given by 


U* 


_ IyJat 
M- 


where AT=T max -T m i n represents the maximum temperature difference at the surface, lyl is the 
absolute value of the derivative of the surface tension with respect to temperature, and 4 is the 
dynamic viscosity. 


5 


g 


t 



Figure 2. The Model Floating Zone 

The temperature maximum in a floating zone occurs between the two ends of the zone. 
We shall take T max to be Too(0), and T m j n to be T 0 , where A = L 0 is the aspect ratio. 

The non-dimensional pressure is 


P = 


p* + ppg*(t) z 
PoU* 2 


Ro 


where p* is the dimensional pressure, g*(t) is the residual gravitational acceleration, z is the 
dimensionless axial coordinate and po is the density corresponding to the reference temperature. 
The temperature is rendered dimensionless using T max - T m j n . With these scales the 
dimensionless equations in a cylindrical coordinate system can be written as 
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I 

r 


d(ru) 

dr 



= 0 , 


( 1 ) 


du du du dp i /d 2 u i du d 2 u u \ 

dt dr dz dr Re \ d r 2 r d r ^ z 2 r 2 f 


(2) 



dw dw dp i /d 2 w | dw d 2 w 

dr dz dz Re \ 0 r 2 r dr d z 2 


+ -GLTg(t), 
Re 2 


(3) 


dT dT dT 

-=T- + U ^r-+ W 

dt dr dz 


j_|^T +i dT + OT| 

M a \dr 2 r dr 0 Z 2 J 


(4) 


where g(t) is a time-dependent dimensionless residual acceleration, and the Reynolds number. 
Re, Marangoni number. Ma and Grashof number Gr are, respectively, 


Re = 


RqU* 

V 


Ma = 


YlATRo 

[IK 


Cr |g*|PATR3 

v 2 


Here, v is the kinematic viscosity, k is the thermal diffusivity, p is the volume thermal expansion 
and coefficient, and |g*| is the (maximum) residual acceleration magnitude 

Boundary Conditions 

At the ends of the zone, the boundary conditions are 

u = T = 0, w = w + (t), at z = + y, (5) 


where w + (t) is zero if the two rigid disks vibrate in phase (this will admit only odd mode 
deformations of the zone surface [ 1 ] ). The symmetry conditions at the centerline r = 0 are 


u = 


dw 

dr 



( 6 ) 


The boundary conditions at the free surface r=R(z) take the form 
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Where 


p- Gg(t)z + X = 


2R~e 


1+ 


3R) 2 
dz 


du /dR\ 2 dw dR /dw du\ 
dr \ dz) dz dz \ dr + dz) 


Re' 1 (C ( V 

-t) 

i i 

(71V 

di) 

a : R 

L /dR \ 2 

1-h-=H 

\dz) 

3/2 

L 

R 

dz 2 


m 


/dR' 7 

I dz 


Id u dw 

fe + lT l+ 


? 3R/3u dw\_ , /dR\ 2 

dzjdr'dzl '■ 3dzJ 


2 /dT dR dT 
dz . 1 dz dr 


( 8 ) 


dR dR „ 

+ U + W TT— = 0 . 

dt dz 


( 9 ) 


1 1 

(dT 

I 1 + 

Or I 

1 dz 

f 

i 

i! 

|2 

\dr 


( 10 ) 


p;-hRo 


C 0 = - , Bi = ^, G = 

To K 


_ to 


Ro 


u * 2 


and are, respectively, the capillary number, Biot number, dimensionless gravitational 
acceleration and y 0 is the mean surface tension and k is thermal conductivity. The force balance 
at the free surface in the normal and tangential directions are given by eqs. (7) and (8), 
respectively. Equation (9) is the kinematic boundary condition at the liquid-gas interface. The 
thermal boundary condition at the interface is given by equation (10) in which the equivalent 
heat transfer coefficient h contains the effect of the radiant and convective heat transfer. The 
constant X in (7) represents a dimensionless reference pressure difference 2 across the interface. 
In model liquid bridge systems with fixed rigid endwalls such as the one discussed here X is 
determined bv the following constant volume constraint 


V = I tc R 2 (z)dz = V 0 = constant . (11) 


Finally, the condition that the contact lines between the liquid end disks are fixed is 
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R = 1 at z = ± A/2. (12) 

Numerical method 

The governing equations are recast in terms of a vorticity stream function together with a 
non-orthogonal coordinate transformation. The latter allows an irregular free "boundary to 
coincide with a coordinate line (or surface) without the need to solve a coupled set of Laplace 
equations. The resulting equations are discretized following a pseudo-unsteady semi-implicit 
difference scheme and solved by the A.D.I. method. The combination of the above methods 
provides a very accurate and economical solution procedure. Generally, the tangential stress 
condition and kinematic condition at the free surface are solved along with the Navier-Stokes and 
continuity equations. The normal force balance condition is used to determine the free surface 
shape. In addition, an "outer" iterative procedure is needed to locate the free surface. 

The unsteady free boundary problem for cylindrical liquid zone is solved as follows. First 
the initial conditions are specified. We shall use initial conditions corresponding to either zero 
residual acceleration or a steady residual acceleration. These states will be calculate using a 
method similar to that described below but designed specifically for the computation of steady 
solutions. A Picard iterative procedure [8] is then followed: 

1. Guess the free surface shape for the new time step; 

2. Obtain the approximate temperature and velocity fields by transforming the governing 
equations and boundary conditions to a circular cylindrical domain via a non-orthogonal 
transformation and solve them using a pseudo-unsteady semi-implicit method; 

3. Obtain the pressure at the free surface by integrating the transformed momentum equation; 

4. Use the normal force balance condition at free surface to decide how to update the free 
surface location; 

5. Return to step 2. Repeat until convergence is obtained by satisfying all equations and 
boundary conditions to a specified degree of accuracy for this time step. 


Non-orthogonal Transformation 

The region occupied by the liquid zone is transformed into a fixed cylindrical 
computational region (see Fig. 2) using a non-orthogonal coordinate transformation, i.e. 


n = z , 



(17) 


This transforms the domain 


- < z < 0 < r < R(z), 

onto the rectangle. 

-A< Z <A 0< r < R(z), 

It then follows that 


9 


3 ? i a a ) a _ c 3 r a 

9r R ^ 9z 9r| R 9r| ^ 


The advantage of this transformation is that the free boundary coincides exactly with a 
coordinate line in a computational grid and regeneration of mesh during the outer iteration is 
avoided. 

We have checked our method by comparing our results with those of Hyer et al. [9] has 
been made. The convergence of the solution was also checked by varying the spatial resolution 
of the mesh. This was particularly important at higher values of Re, where, owing to the space- 
centered-differences [10] the stream function was prone to exhibit “wiggles” if the grid Peclet 
number exceeded 2 in the vicinity of the disks. As observed by Ryskin and Leal [11], an 
increase in the mesh resolution was found to eliminate this problem. Thus, the use of upwind 
differencing, and consequent loss of accuracy [ 10,1 1], to avoid this problem was not necessary. 

4. Results 

The results of our calculations are summarized in this section. More details and the 
results of our analysis of the HISA data are given in the appendices. 

The sensitivity of an indium float zone was analyzed using the 2D model described 
earlier. The input accelerations used ranged from a simple sinusoidal disturbance to a disturbance 
synthesized from an the Honeywell In-space accelerometer. Our studies indicate that the 
components of the disturbance with frequencies around 2.5 -4 Hz results in the greatest response 
for any given acceleration amplitude. This is because the eigenfrequency of the zone examined is 
3.4 Hz. Figure 6 in appendix I shows a set of sensitivity curves for an indium zone which shows 
the enhanced sensitivity of the bridge to vibration in the vicinity of the bridge eigenfrequencies. 
This result is in qualitative agreement with surface motions observed on video taken of the zone 
on STS-52, and also agrees with our results simulated from HISA data which shows a peak 
acceleration associated with 3 frequencies in the 1-5 Hz range. The largest peak corresponds to 
2.7 Hz and is associated with use of the treadmill (see appendix II). The other two peaks are 
associated with structural modes and correspond to 1. 8 Hz and 3.5 Hz. The maximum tolerable 
amplitude to produce a significant observable deformation of the bridge surface was exceeded 
about 39% of the time that HISA was operational. More details on the acceleration history are 
given in appendix II. 


5. Publications, papers presented at conferences and students graduated 

Publications 

J.I.D. Alexander, H. Cordier, J. Ouazzani and Y. Zhang, "Numerical Simulation of 
Thermocapillary flow under Zero and Low Gravity Conditions," AIAA Paper no. 93- 
0254. 31st Aerospace Sciences Meeting. Reno. Nevada. 1993. 

M.J.B. Rogers, J.I.D. Alexander and J. Schoess, Detailed Analysis of the Honeywell In-space 
Accelerometer- STS 52," Microgravitv Science and Technology VI. 28-33. 1993. 

J.I.D. Alexander and Y.Q. Zhang, “The Sensitivity of a Nonisothermal Liquid Bridge to Residual 
Acceleration,” Microgravitv Science and Technology . IV, 128-129, 1991. *extended 
abstract; Full article in: Microgravitv Fluid Mechanics , ed. H.J. Rath, Proceedings of the 
IUTAM Symposium, Bremen, Germany (Springer Verlag, Berlin, 1991) pp. 167-174. 

Papers presented 

“Analysis of Experiment Sensitivity to Residual Acceleration”, International Symposium on 
Microgravitv Science and Application . Beijing, China, May 10-12, 1993 
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“Modelling or Muddling? Analysis of Buoyancy Effects on Transport under Low Gravity 
Conditions,” presented at the World Space Congress . Washington D.C. August 28- 
September 5, 1992. 

“Numerical Simulation of Thermocapillary Flow Under Zero and Low Gravity Conditions” 
presented at the AIAA 31st Annual Aerospace Science Meeting . Reno, Nevada, USA. 
Jan. 11-14, 1993 

“Detailed Analysis of the Honeywell In-Space Accelerometer Data STS-32”, (poster) presented 
at the World Space Congress . Washington D.C. August 28-September 5, 1992. 

“The Sensitivity of a Nonisothermal Liquid Bridge to Residual Acceleration," 5th Annual 
Alabama Materials Research Conference . September 25-26 1991, Birmingham, Alabama. 
"The Sensitivity of a Nonisothermal Liquid Bridge to Residual Acceleration," IUTAM 
Symposium on Microgravity Fluid Mechanics . September 2-6, 1991 Bremen. Germany. 
“Numerical Analysis of the Effects of Low Gravity on Convection and Transport," Case Western 
Reserve University, Department of Mechanical Engineering . February 6th. 1999 
"Numerical Simulation of Low-g Fluid Transport", AIAA Short Course on Low-Gravity 
Fluid Mechanics. January 10-12. 1992. Reno. Nevada ~ 

Students graduated 

Helene Cordier, M.S. Mechanical Engineering, 1993. 
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The sensitivity of a non-isothermal liquid bridge to residual 
acceleration 


J. Iwan D. Alexander and Yiqiang Zhang 

Center for Microgravity and Materials Research 

University of Alabama in Huntsville, Huntsville, Alabama 35899, USA 

Abstract 

Liquid bridges appear in a variety of industrial processes, for example in the well-known floating-zone 
crystal growth technique. This crystal growth method has received much attention in recent years. 
In particular, there have been a variety of experiments on spacelab missions. These experiments are 
motivated by the fact that the microgravity environment affords the possibility of an increase in the 
stability of the melt meniscus and a reduction in buoyancy-driven convection. However, within the 
spacecraft there is a residual acceleration with variable magnitude and orientation. Under certain 
conditions, the response of the free surface of a liquid bridge to time-dependent residual accelerations 
will lead to zone breakage. In this paper the steady and unsteady behavior of isothermal and non- 
isothermal liquid bridge systems under normal and low gravity conditions is examined. The full non- 
linear governing equations are recast in terms of a stream-function vorticity formulation together with 
a non-orthogonal coordinate transformation. The latter allows an irregular free boundary to coincide 
with a coordinate line (or surface) without the need to solve a coupled set of Laplace equations. The 
resulting equations are discretized using a centered finite difference scheme for space, and an Adams- 
Bashforth-Crank-Nicolson scheme is used for time. The equations are solved by the A.D.I. method 
and a Picard type iteration is used on the boundary condition for the balance of force normal to the 
free surface. For non-isothermal bridges, residual acceleration affects the system by causing internal 
buoyancy flows and fluctuations in the shape of the bridge which interact with the thermocapillary 
flow caused by surface tension gradients. For the cases examined, the shape of the bridge is found 
to be more sensitive to typical spacecraft accelerations than the buoyancy driven flow. The effect of 
thermocapillary flow on the surface shape is found to be small for the range of capillary and Reynolds 

1. Introduction 

Liquid bridges appear in a variety of industrial processes, for example the well-known 
floating- zone crystal growth technique [1]. This crystal growth method has received much attention 
in recent years [2-10]. In particular, there have been several related experiments on spacelab missions 
[ 1 1-14]. These experiments are motivated by the fact that the microgravity environment affords the 
possibility of an increase in the stability of the melt meniscus and a reduction in buoyancy-driven 
convection. However, within the spacecraft there is a residual acceleration with variable magnitude 
and orientation. Under certain conditions, the response of the free surface of an isothermal liquid brid ge 
to time-dependent residual accelerations will lead to zone breakage [ 1 3- 1 5] . In this paper we examine 
the interaction between convection caused by the response of the free surface of the zone to oscillatory 
axial residual acceleration and convection due to thermocapillary and internal buoyancy forces. 



L 


Fig. 



Gr = 0.98, A = 4 and Pr = 0.0127. 


The model liquid bridge. 


2. Formulation 

2.1 Description of the Model 

A cylindrical liquid zone (see Fig. 1 ) is contained between two parallel coaxial circular rigid 
disks (radius = Ro) separated by adistance L. The liquid is a non-isothermal Newtonian fluid for which 
the Boussinesq approximation holds. The bridge is held between the disks by surface tension. The 
free surface of the bridge is a gas-liquid interface and is described by r =R(z,t). Each disk is maintained 
at a constant temperature To. Surface heating is provided through a parabolic function, Too, which is 
a function of the axial coordinate and is an approximation to the heating profile associated with typical 
floating zone crystal growth experiments. The heat transfer coefficient at the free surface is denoted 
by h. In addition, we make the assumptions that the residual acceleration is parallel to the cylinder 
axis, the velocity and temperature field and the deformation of free surface are axisymmetric, and we 
take the surface tension at the free surface to be a linear function of the temperature. Motion of the 
end disks perpendicular to the axis may also occur in practice. The restriction to axial acceleration 
precludes an analysis of the effects of such motions. 

The governing equations are made dimensionless by scaling length, time and velocity with 
Ro. Ro/E and U*, respectively. Here U* is a characteristic velocity given by 



where AT = T MAX - T M[N represents the maximum temperature difference along the surface, kj is 
the absolute value of the derivative of the surface tension with respect to temperature, and ji is the 



dynamic viscosity. The difference AT is used to non-dimensionalize temperature. 

The temperature maximum in a floating zone occurs between the two ends of the zone. We 
shall take the dimensionless ambient temperature Tmax to beToo(O), and Tmin to beT«,(±A/2), where 
A = L/Rq is the aspect ratio. 

The non-dimensional pressure is 


P* +pog*(t) z 

p= iq — Ro* 


PoU 

where p* is the dimensional pressure, z is the dimensionless axial coordinate, po is the density cor- 


responding to the reference temperature and g*(t) = go + gisin(2/t/t) is the residual gravitational 
acceleration. 


2.2 Basic Equations 

With the scales presented in 2.1 the governing dimensionless equations in a cylindrical 
coordinate system can be written as 


1 9(ru) 9w 
r 9r 9z 


= 0, 


( 1 ) 


9u 9u 9u 9p l / i 9u 9^ u ] 

9t 9r 9z 9r Re \ ^ r 9r 9z 2 r 2 / 


9w 9w dw dp 1 / 9 i 9w 

9t 9r 9z 9z Re \ ^ r 9r 9z 2 


+-OLTg(t), 

Re 2 


(3) 


9T 9T 9T i / 9^ i 9T 9 t\ 

— + u— — + w— — = — +— — — i , 

9t 9r 9z Ma \ ^ r 9r 0 Z 2 J 

where g(t) = g*(t)/go is a time-dependent dimensionless residual acceleration, and 


(4) 


Rt = M± , Ma Ji& Gr J«, 

V |IK V 2 

are, respectively, the Reynolds number, Marangoni number and Grashof number. Here, v is the ki- 
nematic viscosity, k is the thermal diffusivity, (3 is the volume thermal expansion and coefficient. 
The boundary conditions at the rigid end disks are 


u= T = 0, w = w 1 (t), at z = ± 


(5) 


where w ± (t) is zero if the two disks vibrate in phase (this will admit only odd mode deformations of 
the zone surface [15] ). The symmetry conditions at the centerline r = 0 are 

9w 9T a 

u- aF-*" a 

The boundary conditions at the free surface r = R(z) take the form 


(6) 





Fig. 3. The instantaneous dimensionless stream-functions with Re = 2899, Gr = 0.98, A = 4 
and Pr = 0.0127, at (a) 0.52 s, (b) 0.59 s, (c) 1 s, (d) 1.5 s (e) 1 .56 s, (f) 2 s, after application of 
an additional 2.5x10“ 2 g, 0.5 Hz, axial acceleration. 
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where the capillary number, Biot number and Froude number are, respectively. 
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and y 0 is the mean surface tension and k is the thermal conductivity. The force balance conditions 
normal and tangent to the free surface are given by eqs. (7) and (8) respectively. Equation (9) is the 
kinematic boundary condition at the liquid-gas interface. The thermal boundary condition at the 
interface is given by equation ( 1 0) in which theequivalent heat transfer coefficient, h, models the effect 
of the radiant and convective heat transfer between the bridge and the surrounding environment. The 
constant X in (7) represents a dimensionless reference pressure difference across the interface which 
for this system is determined by the following constant volume constraint [9,16] 

rf 

V = I nR^zjdz =V 0 = constant . (11) 

J=A 

2 

Finally, the condition that the contact lines between the liquid end disks are fixed is 

R = 1 at z = ± A/2. ( 1 2) 



Fig. 4. The dimensionless temperature field with Re = 2899, Gr = 0.98, A = 4 and 
Pr = 0.0127, at (a) 0.52 s, (b) 0.59 s, (c) 1 s, (d) 1.5 s (e) 1.56 s, (0 2 s, after 
application of an additional 2.5x1 0"^ g, 0.5 Hz, axial acceleration. 




2.3 Solution Method 

In the present investigation, the governing equations are recast in terms of a stream-function 
vorticity formulation. The stream-function is defined by 


u = - 


dz 


1<W 

r dr ' 


(13) 


A non-orthogonal coordinate transformation, 

ri= z 


4 = 


(14) 


R(z,t) * 

allows an irregular free boundary to coincide with a cylindrical coordinate line (or surface) without the 
need to solve a coupled set of Laplace equations [17,18]. The resulting equations are discretized 
following a semi-implicit difference scheme and solved by the A.D.I. method. The conditions for force 
balance tangent to the surface and kinematic condition at the free surface are solved along with the 
Navier-Stokes and continuity equations. The condition for the force balance normal to the surface is 
used together with an “outer” iterative procedure to determine the free surface shape. 

The unsteady free boundary problem for a cylindrical liquid zone is solved as follows. The 




S/Vmin = "3.078x 1 0' 
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Fig. 5. The instantaneous dimensionless stream- 
function, with Re = 0, Gr = 0, A = 4 and Pr = 0.0127, 
at (a) 0.59 s, (b) 1 s, (d) 1 .5 s (e) 2 s, after application 
of an additional 2. 5x10' 2 g, 0.5 Hz, axial acceleration. 
Note that for this ^ase only, v P= ‘P/kRo , where 1 ? is the 
dimensional stream-function cm 3 s 1 . 



Fig. 6. Tolerable acceleration [g] vs. Frequency [Hz] 
predicted by the ID model [15] for (1) breakage and 
(2) 10% shape change, for an isothermal liquid zone 
with A = 4 and the thermophysical properties of 
indium. 




initial conditions correspond to either zero or finite steady residual acceleration situations with a steady 
thermocapillary flow. These states are calculated using a method close to that described below [16]. 
For the unsteady calculation our solution scheme is similar to that used by Kang and Leal [17] and 
Ryskin and Leal [18]. The following Picard iterative procedure [19] is adopted: 

1. guess the free surface shape for the new timestep; 

2. obtain the approximate temperature and stream-function, vorticity and velocity fields by transform- 
ing the governing equations and boundary conditions to a circular cylindrical domain via a non- 
orthogonal transformation and solve them using a semi-implicit method; 

3. obtain the pressure at the free surface by integrating the transformed momentum equation; 

4. use the condition for the balance of force normal to the free surface to decide how to update the free 
surface location; 

5. return to step 2. Repeat until convergence is obtained by satisfying all equations and boundary 
conditions to a specified degree of accuracy for this timestep. 

3. Results and Discussion 

The following results were obtained for a liquid zone corresponding to the physical properties 
of molten indium subject to an axial acceleration with a frequency of 0.5 Hz. Fig. 2 depicts the initial 
state of the system . A steady axial acceleration of magnitude 10" 4 g ( 10* 3 m s' 2 ) acts along the negative 
z-direction. Two equidimensional toroidal rolls indicate that surface-driven flow is dominant (Re = 
2899, Gr = 0.98). The isotherm distribution shows that heat transfer is mainly by conduction, although 
some distortion of the isotherms by the flow is evident. For indium, Lind [20] has reported that the 
surface tension increases with increasing temperature i.e. Yr>0- (Note that surface contamination may 
have affected the measured temperature dependence of surface tension in this case.) Thus, the flow 
direction at the surface is toward the center (i.e. the higher temperature region). Figs. 3 and 4 illustrate 
the etfect of an additional acceleration component which varies sinusoidally with a frequency of 0.5 
Hz. Figure 5 depicts the response of the zone to the same disturbance, but with Gr = Ma = 0. (Note 
that, for this case, k/L rather than AT/p was used as the velocity scale.) In another case with Re = 
2899, Gr = 0.98 and the surface constrained to be a circular cylinder, no observable response occurred. 
Clearly the system is more sensitive to the effects of free surface motion than internal buoyancy. Given 
that the value of steady acceleration used is extreme for spacecraft acceleration environments [2 1 ] , we 
may conclude that for systems where internal buoyancy-driven effects are not manifested (in this case 
because they are swamped by the surface-driven flow) it suffices to examine the response of the free 
surface only. Furthermore, a comparison of our full axisymmetric results with those obtained with a 
simplified ID isothermal model indicates that, at least for the conditions examined, the ID model may 
be used to reliably predict liquid zone (isothermal and nonisothermal) sensitivity. The ID model is 
described in detail in [15]. Fig. 6 shows the sensitivity of an indium liquid zone to axial acceleration. 
The curves are based on results obtained using two sensitivity criteria. The first is determined by 



breakage of the bridge, the second is whenever the bridge shape changes by more than 1 0% of its static 

shape, i.e. R(z,t) - R(z,0) = .lR(z,0). 
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Detailed Analysis of Honeywell In-space 
Accelerometer Data - STS-32 ***■ 

<? 3 // 3^070 


The Honeywell In -Space Accelerometer system (HIS A ) col- 
lected data in the mid-deck area of the Shuttle Columbia 
during the flight of STS-32. January 1990. The resulting data 
were to he used to investigate the response of crystal mi- 
crostructure to different types of residual acceleration. The 
HISA , using Sundstrand Q A -2000 sensors . is designed to 
detect and record transient and oscillatory accelerations. The 
sampling and electronics package stored averaged accelera- 
tions over two sampling periods; two sampling rates were 
available: I Hz and 50 Hz. Analysis of the HISA data fol- 
lowed the CMMR Acceleration Data Processing Guide, con- 
sidering in -house computer modelling of a float -zone indium 
crystal growth experiment. Characteristic examples of HISA 
data showing the response to the primary reaction control 
system (PRC'S). Orbiter Maneuvering System (OMS) opera- 
tions and crew treadmill activity are presented. Various ar- 
biter structural modes are excited by these and other 
activities. Such modes are in the l Hz to 20 Hz range. Of 
particular note is the distinct frequency of crew footfalls 
during exercise periods and the resultant excitation of arbiter 
frequencies. Acceleration vector magnitudes ranged from 
10 4 g to >10 - g. Acceleration values recorded at the crys- 
tal growth site exceeded time domain tolerance limits in 
approximately 35 % of the data . especially during PRCS 
firings , OMS burns , and treadmill exercise. Frequency do- 
main limits , however* were not exceeded. 

J Introduction 

The Honeywell /n-Space Accelerometer system (HISA) 
collected data on the mid-deck of the Space Shuttle Colum- 
bia during the flight of STS-32. January 1990. The HISA 
was part of the A/icrogravity Disturbances Experiment 
(MDE) which also included the Rockwell Fluids Experi- 
ment Apparatus (FEA). The FEA is a float-zone materials 
processing furnace; for this flight, it was used to grow single 
crystals of metallic indium. At various times during growth, 
acceleration sources were activated in order to investigate 
the response of crystal microstructure to different types of 
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residual acceleration. The HISA acceleration measurements 
are presented and discussed herein. 

In sect. 2 we outline the MDE location and HISA 
operation in the orbiter mid-deck. Sect. 3 discusses our 
processing plan, based on estimated acceleration sensitivi- 
ties. In sect. 4 we present and discuss time history and 
frequency domain plots of segments of the HISA data. 
These data are compared to data from other missions in 
sect. 5 and suggestions are made for revisions of the HISA 
system before its further use in low-gravity environments. 

2 Experiment Set-up 

We present here a brief overview of the MDE; more de- 
tailed discussions can be found in Schoess [1] and Dunbar et 
al. [2], During operation the MDE was mounted at locker 
MA9F in the at t mid-deck of Columbia. Crew members 
observed and videotaped the shape changes of the liquid 
indium zone. Post flight crystal characterization was 
planned to investigate crystallographic defects and impurity 
content; the experiment was designed to investigate the 
relationships between crystal quality and variations in ac- 
celeration level, growth rate, and zone stability [1]. 

The HISA, composed of three Sundstrand Q A -2000 sen- 
sors and an electronics package, is designed to detect and 
record transient and oscillatory accelerations. Table 1 lists 
HISA performance specifications. Data were recorded at 
both 1 Hz and 50 Hz sampling frequencies as averages of 
acceleration impulses over two sampling periods. They do 
not, therefore, represent instantaneous acceleration. The 
electronics package is such that positive and negative aver- 


Table 1. HISA performance specifications , after Schoess [I] 


parameter 

performance 

orientation 

three-axis orthogonal 

range 

10 A 10 10 'g at 1 Hz 
10- 5 to 10 at 50 Hz 

accuracy 

resolution 

± ( 1 % |reading| + 0.00002) ? 
< 1.0 x 10-\ ? at 1 Hz 


8.7 x 10 - '" g at 50 Hz 

frequency response ( ±5%) 

0.025 to 19.500 Hz 

DC bias 

none (AC output) 

sample data rate 

50 Hz. 1 Hz 

communications 

RS-422/ASCII format 

size 

20.3 x 9.7 x 5.3 cm ( 1 .046 cm 5 ) 

mass 

1.8 kg 

power 

5.6 W (28 V) 
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ages were recorded separately. This recording technique 
resulted in twelve data files for each time period studied: 
±.v, ± v, and ±z data for both the 1 Hz and 50 Hz sam- 
pling rates. 

3 Data Analysis 

The HISA data files were provided by the NASA .4ccerler- 
ation Characterization and Analysis Project (ACAP). The 
data were in tabular form with the averaged acceleration 
values and sensor temperatures listed for each time mark. 
The temperature variation recorded by the system was not 
enough to warrant temperature corrections to the data [3]. 
A problem with the recording hardware responsible for the 
switch between the two sampling rates caused spurious 
values to be recorded approximately once a second in the 
50 Hz data. We tested three methods to resolve this prob- 
lem: replacing the spurious points with zero, with a “not-a- 
number" marker, and with the same value as the preceding 
data point. The three techniques showed similar results, and 
we chose to use the zero replacement method. 

Our approach to the analysis of the HISA data follows 
the Acceleration Data Processing Guide developed at the 
Center for Microgravity and Materials Research. This 
guide was developed to assist low-g experimenters inter- 
ested in accelerometer data in determining their data access 
and processing needs based on their knowledge of experi- 
ment and acceleration sensor location and of experiment 
sensitivity to acceleration [4]. We assessed the sensitivity of 
the float-zone indium crystal growth experiment based on 
computer modelling of the experiment [5 -6]. Table 2 lists 
the sensitivity limits used in our analysis. 


Table 2. Sensitivity limits used in analysis, after Alexander and 
Zhan% [6] 


frequency and magnitude ranges 

10“ J £ for 10~ 2 to 0.1 Hz 1 

of interest 



10 ' (J for 0.1 to 10 Hz 

maximum tolerable acceleration 

io-\? 

{ steady) 


maximum tolerable acceleration 


(transient) j 

10-\? | 


Typically, sensitivity limits are used to determine what 
segments of an acceleration data base to analyze in con- 
junction with experimental results. Sections of accelerome- 
ter data in which the sensitivity limits are exceeded are 
selected using a threshold detection routine. A relatively 
limited amount of accelerometer data was recorded by the 
HISA during STS-32 (^10 megabytes representing ^5 
discontinuous hours). Because of this, the entire data base 
was analyzed, with increased attention given to segments 
where the sensitivity limits were exceeded and where known 
acceleration sources were active. 

A typical processing routine for residual acceleration 
data uses the basic techniques of time series and spectral 
analysis. These include calculating series statistics (mean, 
minimum, maximum, variance), plotting time histories of 
individual axes of data and vector magnitudes, and evaluat- 
ing and plotting amplitude and power spectral densities to 


investigate the frequency components present and their 
relative strengths [4, 7-9]. 

As mentioned in sect. 2, the HISA data were digitized so 
that positive and negative averages were recorded sepa- 
rately. This recording technique necessitated some special- 
ized analysis. Because we are generally more interested in 
the magnitude of the acceleration vector at a location than 
in the individual axial acceleration readings, we used the six 
data files from each sampling frequency to form the vector 
magnitude, that is |a| = x / ( ±a 2 x + ±a 2 v -f ±al), where 
a = ( ± ci x , ±a v , ±a : ) represents the positive and negative 
averages of the accleration vector. 

We transformed the HISA time histories into the fre- 
quency domain in a similar manner. Each of the six data 
files for each sampling rate was Fourier transformed to 
create the corresponding amplitude spectrum. To demon- 
strate that this technique is appropriate, we made a visual 
comparison of the spectra with spectra resulting from more 
standard data bases. For a more quantitative test, we 
processed data from the STS-40 flight of the Space Acce ler- 
ation Measurement System (SAMS) to simulate the HISA 
sampling and averaging technique. The time and frequency 
domain magnitudes obtained with the resampled SAMS 
data were comparable to those obtained through our stan- 
dard SAMS processing techniques. 

The visual comparisons of the HISA results to other 
data bases and the SAMS simulated data show that the 
typical structural modes seen in orbiter accelerometer data 
were present in the HISA data. The lowest frequencies, 
unfortunately, are obscured to some extent by a remnant 
effect of the “ + " digitization. The lowest frequency compo- 
nent resulting from a Fourier transformation of a time 
series is a representation of the mean value of the series. 
The ** -f- digitization of the HISA data creates an artificial 
non-zero mean that corrupts the mean represented in 
the corresponding spectrum. Because of this effect, we gen- 
erally disregard the lowest spectral components calculated. 
For the shortest HISA data liles (the various thruster 
firings), this results in a lower frequency limit of about 
0.1 Hz. 

4 HfSA Data Plots 

We present several characteristic examples from the HISA 
data base. Fig. 1 shows a time history and amplitude 
spectrum for 50 Hz data collected during a primary reac- 
tion control system firing. A high thrust pulse occurs about 

5 s into the record, with possible lower magnitude thrusts 
occuring at about 45 and 90 s. The amplitude spectrum 
indicates that several orbiter structural modes are excited 
during this period, most notably in the 3, 5, 7, 10, and 
13 Hz regions. The 17 Hz component, which is related to 
the Ku band antenna dither and a nearby orbiter structural 
mode, is excited. The lower magnitude, narrow band signals 
at 16 and 20 Hz are unusual, and at present their origin is 
uncertain. Early analysis of SAMS mid-deck data from 
STS-43, August 1991, suggests that the 20 Hz frequency 
may be a local, mid-deck structural or equipment-related 
excitation. Fig. 2 shows comparable plots for the 1 Hz data. 
Note that the amplitude spectrum plot can be used to 
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Fig. I. Example of a primary reaction control system firing recorded hy the HIS A in the mid-deck of Columbia on STS-32 , 50 Hz data , 
(a) acceleration rector magnitude , (b) combined amplitude spectrum 



supplement the lower frequency information obscured in 
fig. 1. Data recorded during an Orbiter A/aneuvering Sys- 
tem (OMS) burn are shown in fig. 3. The burn produces a 
short pulse-like acceleration response at the mid-deck loca- 
tion of the HISA. Structural modes in the 3.5 and 7 Hz 
range are excited. 

Figs. 4-6 show HISA data collected during a period of 
crew exercise on the mid-deck treadmill. The treadmill was 
located approximately 2 m from the HISA. Fig. 4 shows a 
1 Hz data window, and corresponding amplitude spectrum, 
collected during a crew treadmill exercise period. Note the 
high and low magnitude segments related to different exer- 
cise intensities. Fig. 5 shows 50 Hz data collected during a 
running or jogging (high intensity) exercise phase. Fig. 6 
shows a segment from a less intense exercise stage, such as 
a walking or cool down period. Note that the overall 
acceleration level of fig. 6 is lower in magnitude than that of 
fig. 5. The most interesting aspect of these plots is the 
difference in the amplitude spectra. The spectra alone could 
be used to identify that the crew member was exercising at 
different levels of intensity. 


This is explained as follows. The higher magnitude (run- 
ning) segments of the treadmill data are dominated by a 
2.7 Hz component, fig. 5. This frequecy is related to an 
orbiter structural mode and is being driven by the footfalls 
of the crew member. In fig. 6, a slower exercise rate is 
represented (walking) and has a correspondingly lower 
(requency of footfalls. The 1.8 Hz component may be a 
higher harmonic of the 0.9 Hz component, or the two 
components may represent the crew member's upward and 
downward movements which differ in magnitude due to the 
treadmill restraint system. In this case, the 3.5 Hz structural 
mode is excited, possibly related to a doubling of the 1.8 Hz 
mode. 

Because of the relatively small amount of data collected 
with the HISA during STS-32, we analyzed the entire data 
set, rather than limit the analysis based on sensitivity values 
given in table 2. Comparison of the HISA time histories to 
the limits in table 2 indicates that the transient time domain 
limit of 10 was exceeded about 35% of the time as 
recorded at 50 Hz and 39 % of the time at 1 Hz. The 
transient limit typically was exceeded during periods of 
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fig J Example of an arbiter maneuvering system firing recorded by the HISA in the mid-deck of Columbia on STS-32. 50 Hz data . 
ft;) acceleration magnitude , f/>) combined amplitude spectrum 



Fig. 4. Example of acceleration level during crew treadmill exercise period recorded by the HISA in the mid-deck of Columbia on STS-32 , 
/ //r data. Note the different acceleration levels corresponding to different stages of the treadmill exercise protocol 



(a) TIME [sec) (b) FREQUENCY [Hz] 

Fig. 5. Segment of 50 Hz data from higher magnitude section of fig. 3 , (a) acceleration vector magnitude . (h) combined amplitude spectrum. 
Note the 2.7 Hz structural mode excited by footfall frequency 
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Fig. 6 Segment of 50 Hz lower magnitude section of fig. 5. (a ) acceleration vector magnitude, (h) combined amplitude spectrum. Note the lower 
footfall fretfuencv 


crew exercise and orbiter maneuvers. The sampling frequen- 
cies used, the lengths of the time windows recorded, the 
limited amount of data storage capacity, and the analysis 
limitations related to the digitization precluded the 

testing of the steady-state time domain limit. Note that the 
HISA was designed only to measure if-jitter (transient and 
oscillatory accelerations). 

Evaluation of the frequency domain tolerance limits is 
particularly important in the case of float-zone experiments. 
These experiments involve a molten liquid zone suspended 
between two rods. The free surface of the zone will move if 
subjected to a disturbance. In particular, for vibrations, the 
zones tend to be most sensitive to disturbances with fre- 
quencies corresponding to eigenfrequencies of the bridge. 
For the case of an indium floating zone, these frequencies 
lie in the 0.1 to 10 Hz band [6]. For the HISA data 
analyzed, the maximum tolerable magnitude for this range 
was not exceeded. The 10 ~ 2 to 0.1 Hz frequency range limit 
in table 2 was not exceeded for segments of data during 
treadmill activity and quiet periods. The lowest frequencies 
of this range cannot be evaluated for the segments of 
thruster firings because the short time windows recorded 
impose a lower frequency analysis limit in the 10 _2 Hz 
range. For frequencies above 2-3 • 10 -2 Hz, the sensitivity 
limit for zone breakage was not exceeded. 

5 Discussion 

The results obtained from processing the STS-32 Honeywell 
In-Space Accelerometer data are comparable in magnitude 
and character to the expanding orbiter low-gravity mea- 
surement data base. The HISA was one of the first systems 
to be used in the orbiter mid-deck and. consequently, the 
data should serve as a standard against which later results 
will be compared. 

Acceleration vector magnitudes recorded by the HISA 
for primary reaction control system firings (figs. 1 and 2) 
are consistent with data recorded in the Spacelab module 
during the Spacelab 3 and Spacelab Life Sciences 1 missions 
and in the payload bay during STS-32 [7, 11, 12]. It should 
be noted, however, that in some cases the accelerations 


exceed the accelerometer maximum range, so we do not 
know the actual maximum level of acceleration related to 
thruster firings. 

As can be seen in all these data bases, the frequency 
components of a time segment of acceleration data are 
generally dominated by the various orbiter and Spacelab 
structural modes [4, 7-10]. For the case of acceleration 
data collected during crew treadmill exercise (figs. 4-6), 
vibration frequencies related to the crew member's footfalls 
on the treadmill and subsequent pull on the restraint system 
are not only detectable in spectral representations, but are 
also seen to excite structural modes of the same or related 
frequencies. Analysis of HiRAP data reported in Dunbar et 
al. [11] indicates that the magnitude difference between 
exercise stages is translated to the orbiter payload bay at an 
overall reduced magnitude, but that the footfall frequencies 
are not evident. We interpret this as an indication that only 
the vibration of the excited structural modes propagates the 
distance 18 m between the treadmill and payload bay. The 
footfall source excites structures local to the mid-deck 
treadmill and HISA sensors. 

Despite the fact that the HISA data are consistent with 
other data bases, we think that there are several aspects of 
the system that should be modified before future flights. We 
believe that the “ + " digitization method is a relatively 
undesirable form of data conversion for storage. Modifica- 
tion of the electronics to allow recording of a single time 
series for each of me three orthogonal axes would alleviate 
several processing difficulties. Correction of the sampling 
rate switching problem would also obviate data replace- 
ment or filtering schemes and increase the number of mean- 
ingful data points available to the experimenter. 
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NUMERICAL SIMULATION OF THERMOCAPILLARY FLOW UNDER ZERO 
AND LOW GRAVITY CONDITIONS 


Abstract 

This paper discusses the numerical solution 
methods and results of steady and unsteady thermocapillary 
(surface-tension) and buoyancy driven flows in 2-D cavities 
and liquid columns. The 2D cavity was assumed to be 
square with one free surface with a zero Capillary number 
(i.e. the free surface was constrained to be Hat). A 
pseudospectral method was used to solve steady and unsteady 
surface tension-driven and mixed buoyancy-surface tension 
Hows in a square cavity. For the liquid column a finite- 
difference scheme based on a Picard iteration was used to 
solve for the flow, temperature and free surface shape. The 
surface of the liquid column was allowed to deform and, as 
for the 2D cavity, the surface tension was assumed to depend 
on temperature. 


Nomenclature 

L length of free surface 

N r number of points in the radial direction 

N z number of points in the axial direction 

Pr Prandtl number 

r radial coordinate 

Re surface Reynolds number 

AT temperature difference between the sidewalls of the 
square cavity 
z axial coordinate 

Yj rate of change of surface tension with respect to 
temperature 

p shear viscosity 

v kinematic viscosity 


I. Pseudospectral Method: Square Cavity 
We assumed that a constant viscosity Boussinesq 
fluid was contained in a square cavity with differentially 
heated isothermal sidewalls, a lower adiabatic boundary and a 
flat upper surface at which motion was driven by gradients 
in (temperature dependent) surface tension. A Chebyshev 


pseudo-spectral collocation method^ was used in 
conjunction with an influence matrix technique ^ to solve the 
problem in stream-tunction vorticity form. Results of the 
calculations were carefully compared to results obtained 
previously using finite difference*^ and finite volume 4 
methods. 

For steady flow, results of the comparison are 
presented in Table I for Reynolds numbers (Re=LATy T /pv) 

in the range 1000 < Re<10 4 and for Prandtl numbers. Pr= 
v/k, in the range 10' 2 < Pr < 1. At Re > 1000. 
convergence was readily obtained for a given number of 
collocation points (grid size), but some effort had to be made 
to ensure grid convergence of solutions. This also appears to 
be the case for finite volume and finite difference methods 
applied to the same problem. We found that, up to Re=10 4 , 
41x41 collocation points were necessary to obtain a grid 
converged solution. That is the solution changed by less 
than 1% after increasing the number of collocation points bv 
50%. 

The 2D surtace tension-driven flow was also 
examined for its sensitivity to sinusoidal vibration. For a 
fixed vibration amplitude and frequency, our results show a 
decreased sensitivity with increasing Re. An important 
aspect of this part of our work again concerned the 
importance of using enough spatial resolution to capture the 
physics. We found for the vibration problem that unless a 
sufficiently fine grid was used (in some cases 65x65 
collocation points) spurious oscillations were introduced. 
Spurious oscillations were not found to occur for the same 
vibration imposed on a closed cavity with the same 
temperature gradient. The importance of spatial resolution 
is also seen from our results for thermocapillary flow in 
liquid columns with deformable surfaces. 

IT Liquid Columns: Centered Finite Differences 

We approached the problem of computing steady 
free boundary problem for thermocapillary and buoyancy- 
driven flow in cylindrical liquid bridge using a Picard 
iterative procedure 5 - 6 as since the location of the free liquid 
surafce is a priori unknown: The iterative procedure is as 
follows: 
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Table I: Comparison of the results of our pseudospectral method (AC) with those of Zebib, Homsy an 

Meiburg (ZHM) 3 , Carpenter and Homsy (CH) 4 . 


Re 

AUTHORS 

M xN 

Nu- 

Nu+ 

\j/ x 102 

1000 

ACO 

17x17 

1.915 

mem 



M 

41x41 

1.926 

MM 



?l 

65x65 

1.926 

1.926 



ZHM 

62x54 

1.93 

1.92 

0.479 | 

2000 

ACO 

17x17 

2.481 

2.507 



11 

41x41 

2.474 

2.474 

0.433 


ii 

65x65 

2.474 

2.474 

0.433 


ZHM 

62x54 

2.47 

2.47 

0.424 

3000 

ACO 

17x17 

2.920 

2.921 



n 

41x41 

2.865 

2.866 

| 


ZHM 

62x54 

2.86 

2.85 

0.404 

4000 

ACO 

17x17 

3.288 

3.255 



" 

41x41 

3.172 

3.175 

1 


ZHM 




■ 


ZHM 





5000 

ACO 

17x17 



0.403 


" 

41x41 



0.377 


11 

65x65 

3.435 

3.434 

0.377 


ZHM 

65x65 

■■ 

3.448 

0.370 


it 

65x65 

WBEm 

3.412 

0.366 


ii 

70*60 

3.428 

3.425 

0.369 


it 

80*80 

3.454 

3.447 

0.373 

6000 

ACO 

41*41 | 

3.654 

3.662 

0.363 


ZHM j 

62*54 

3.65 

3.63 

0.350 

7000 

ACO 

17*17 

HM|| 

4.026 



n 

41*41 

| 

3.865 

■h 


n 

65*65 

3.861 

3.862 

I 


ZHM 

62*54 

3.85 

3.83 

■EB9H 

■BVTCBj 

ACO 

17*17 

4.683 




" 

41*41 

4.361 

4.375 



•? 

65*65 

4.360 

4.366 

0.322 


ZHM 

62*54 

4.34 

4.32 

0.305 


CH 

64*64 

4.36 

4.40 

0.323 


2 

















































































1. guess the free surface shape for the initial 

iterate; 

2. obtain the approximate temperature and velocity 
fields by transforming the governing equations and boundary 
conditions to a circular cylindrical (r,z) domain via a non- 
orthogonal transformation and solve them using a pseudo- 
unsteady semi-implicit method; 

3. obtain the pressure at the free surface by 
integrating the transformed momentum equation; 

4. use the normal force balance condition at the 
free surface to decide how to update the free surface location; 

5. return to step 2. Repeat until convergence is 
obtained by satisfying all equations and boundary conditions 
to a specified degree of accuracy. 

The details of this numerical procedure are discussed 
in previous work^. 

Apart from varying the iterative step size, 
convergence of the solution was checked by varying the 
spatial resolution of the mesh. This was particularly 
important at higher values of Re. where, owing to the space- 
centered-differences, the stream function was prone to exhibit 
“wiggles” if the grid Peclet number^ exceeded 2 in the 
vicinity of the disks. The wiggles are present for N r xN z = 
26x61 and are not eliminated until N z = 101. As observed 
by Ryskin and Leal 7 , an increase in the mesh resolution was 
found to eliminate this problem. We attempted to eliminate 
the wiggles and avoid the need for mesh refinement by 
employing second and third order upwind schemes for 
convective terms. We found that while the wiggles were 
certainly eliminated for N r xN z = 26x61. the mesh still 
needed to be refined in order to obtain grid convergence. 
Since the end result was the same, we concluded that the 
centered difference scheme was preferable. A mesh of N r xN z 
= 26x5 1 was found to be sufficient for the results presented 
here with Reynolds numbers in the range 0<Re<10000. 

We have used above method to examine the 
influence of various parameters on momentum and heat 
transport and meniscus shape. In addition, a favorable 
comparison for full and half zones has been made between 
the results obtained using the combined scheme and those of 
Hyer et ah** who employed a finite element scheme. 
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